# Analyses for "A Red Flag for Public Goods? The Correlates of Civil Society Restrictions"
# Date: 01.08.2024
# Author: Hannah Smidt

#### Preamble ####

# Libraries
library(dplyr)
library(lfe)
library(stargazer)
library(RColorBrewer)
library(dotwhisker)


# Clear working space
rm(list=ls())

# Set WD
setwd("C:/Users/HSmidt/Dropbox/PC/Documents/PaperProjects/2023_2_Paper - Restrictions CSOs and global public goods/Replication files")

# Load data
load("Dataset_For_Analysis.RData")

# Delete string country
data$Country <- countrycode::countrycode(data$cown, origin="cown", destination="country.name")
table(data$Country, useNA="always")

#### Figure 1: Show over-time trends in restrictions and DVs ####

png("Fig1.png", res=150)
ggplot(data
       , aes(x=year, y = RESTRICT_COUNT_NOBAN)) +
  geom_smooth(col="red", method="loess") + 
  #geom_point(col="red", size=0.1, position = "jitter") +  
  xlab("Time in years") +
  ylab("Mean count of CSO restrictions")
dev.off()


#### Figures 2-5: Create variable for restriction increase-episodes ####

## Create variable defining episode of increasing restrictions based on V-Dem
mean(na.omit(data$VDEM_CSOrepress_netl5))
sd(na.omit(data$VDEM_CSOrepress_netl5))

data <- data %>% arrange(cown, year) %>%
  group_by(cown) %>%
  mutate(episode_vdem = ifelse(VDEM_CSOrepress_netl5>-0.03827058+0.6119575 & !is.na(VDEM_CSOrepress_netl5),1,0),
         episode_vdem = ifelse(is.na(episode_vdem), 0, episode_vdem),
         
         episode_vdem = ifelse(year<=2020 & lead(episode_vdem,1)==1,1,episode_vdem),
         episode_vdem = ifelse(year<=2020 & lead(episode_vdem,1)==1,1,episode_vdem),
         episode_vdem = ifelse(year<=2020 & lead(episode_vdem,1)==1,1,episode_vdem),
         episode_vdem = ifelse(year<=2020 & lead(episode_vdem,1)==1,1,episode_vdem),
         episode_vdem = ifelse(year<=2020 & lead(episode_vdem,1)==1,1,episode_vdem),
         
         # For graphs by country
         episode_vdem2 = ifelse(VDEM_CSOrepress_l1<VDEM_CSOrepress & 
                                  VDEM_CSOrepress_l2<VDEM_CSOrepress_l1 & 
                                  VDEM_CSOrepress_l3<VDEM_CSOrepress_l2 , 1, 0),
         
         episode_vdem2 = ifelse(year<=2020 & lead(episode_vdem2,1)==1,1,episode_vdem2),
         episode_vdem2 = ifelse(year<=2020 & lead(episode_vdem2,1)==1,1,episode_vdem2),
         episode_vdem2 = ifelse(year<=2020 & lead(episode_vdem2,1)==1,1,episode_vdem2)
  )


## Which countries have increases in restrictions?
unique(data$Country[data$episode_vdem==1])
unique(data$Country[data$episode_vdem2==1])


## Normalize DVs and IVs for graphs
normalize <- function(x){
  x1 <- (x - min(na.omit(x)) ) / ( max(na.omit(x)) - min(na.omit(x)) )
  return(x1)
}

data$VDEM_publicgoodsspend_n <- normalize(data$VDEM_publicgoodsspend)
data$VDEM_corruptionex_n <- normalize(data$VDEM_corruptionex)
data$VDEM_votebuy_n <- normalize(data$VDEM_votebuy)
data$VDEM_progrClientParties_n <- normalize(data$VDEM_progrClientParties)


## Figures for Change in outcome after start of episode of restrictions

data_plot <- subset(data, year>=1990 & year<=2021)
country_episode_start <- data_plot %>% 
  filter(episode_vdem==1) %>%
  group_by(cown) %>% 
  dplyr::summarize(year_min = min(year) )

data_plot2 <- merge(data_plot, country_episode_start, by="cown", all.x=T)
data_plot2$year_after_restr_increase <- data_plot2$year - data_plot2$year_min

# Plots with V-Dem
png("Fig2.png")
ggplot(subset(data_plot2, year_after_restr_increase>=(-10) &
                year_after_restr_increase<=(5))
       , aes(x=year_after_restr_increase, y = VDEM_publicgoodsspend_n)) +
  geom_smooth(col="#6929c4", method="loess") + geom_vline(aes(xintercept = 0), lty = 4) + 
  #geom_point(col="#6929c4", size=0.1, position = "jitter") +  
  xlab("Years from start of episode of mounting restrictions") +
  ylab("Public goods orientation in national spending")
dev.off()

png("Fig3.png")
ggplot(subset(data_plot2, year_after_restr_increase>=(-10) &
                year_after_restr_increase<=(5))
       , aes(x=year_after_restr_increase, y = VDEM_corruptionex_n)) +
  geom_smooth(col="#00539a", method="loess") + geom_vline(aes(xintercept = 0), lty = 4) + 
  xlab("Years from start of episode of mounting restrictions") +
  ylab("Corruption")
dev.off()

png("Fig4.png")
ggplot(subset(data_plot2, year_after_restr_increase>=(-10) &
                year_after_restr_increase<=(5))
       , aes(x=year_after_restr_increase, y = VDEM_votebuy_n)) +
  geom_smooth(col="#005d5d", method="loess") + geom_vline(aes(xintercept = 0), lty = 4) + 
  xlab("Years from start of episode of mounting restrictions") +
  ylab("Vote buying")
dev.off()

png("Fig5.png")
ggplot(subset(data_plot2, year_after_restr_increase>=(-10) &
                year_after_restr_increase<=(5))
       , aes(x=year_after_restr_increase, y = VDEM_progrClientParties_n)) +
  geom_smooth(col="#0e6027", method="loess") + geom_vline(aes(xintercept = 0), lty = 4) + 
  xlab("Years from start of episode of mounting restrictions") +
  ylab("Clientelist parties")
dev.off()





#### Preparation for multivariate analyses ####

# Control variables
controls <- paste0(" + Control_UCDP_ongoing_conflict_netl5 + Control_VDEM_VerticalAccount_netl5 + Control_VDEM_HorizontalAccount_netl5 + Control_WB_controls_GDP_log_netl5 + Control_WB_controls_population_log_netl5")

# Subset data
data <- subset(data, year>=1990 & year<=2021)

# Make sure that restrictions count is NA if data on restrictions is not available or
# if data on restrictions indicate that all CSOs in the country are banned
data$RESTRICT_COUNT_NOBAN_netl5[which( data$year<=2001 | data$year>=2017) ] <- NA
data$RESTRICT_COUNT_NOBAN_netl5[which( data$BANNED==1 ) ] <- NA

# Keep only data without missings on any of the variables used in the analysis
vars_VDEM_UCDP_WB <- c("VDEM_publicgoodsspend", "VDEM_corruptionex", "VDEM_votebuy", "VDEM_progrClientParties",
                       "VDEM_CSOrepress_netl5", "Control_UCDP_ongoing_conflict_netl5", "Control_VDEM_VerticalAccount_netl5", "Control_VDEM_HorizontalAccount_netl5", "Control_WB_controls_GDP_log_netl5", "Control_WB_controls_population_log_netl5")
data <- data[complete.cases(data[vars_VDEM_UCDP_WB]),]

# Create a binary variable for the presence of populist leaders
data$Control_populistleader_l5yrs_bi <- ifelse(data$Control_populistleader_l5yrs>0,1,0)


# Descriptive Stats table
vars <- c("VDEM_publicgoodsspend", "VDEM_corruptionex", "VDEM_votebuy", "VDEM_progrClientParties",
          "Mech_VDEM_engagedSociety_netl5","Mech_shamingINGO_netl5","Mech_CR_protest_netl5_log",
          "RESTRICT_COUNT_NOBAN_netl5","VDEM_CSOrepress_netl5"
          , "Control_VDEM_VerticalAccount_netl5", "Control_VDEM_HorizontalAccount_netl5"
          ,"Control_UCDP_ongoing_conflict_netl5"
          , "Control_WB_controls_GDP_log_netl5", "Control_WB_controls_population_log_netl5"
          ,"Control_populistleader_l5yrs"
          ,"Control_populistleader_l5yrs_bi")

#### Table 1: Descriptive statistics ####

stargazer(as.data.frame(data[,vars]), type="text", out="Table1.txt"
          , title="Summary statistics"
          , covariate.labels = c("DV1: Public goods orientation in national spending"
                                 , "DV2: Executive corruption"
                                 , "DV3: Vote buying"
                                 , "DV4: Clientelist party linkages"
                                 , "M1: Engaged society (sum of net changes, 5yrs)"
                                 , "M2: Intern. shaming (sum, 5yrs logged)"
                                 , "M3: Protest (sum, 5yrs logged)"
                                 , "Restrictions on CSOs (sum of net changes, 5yrs)"
                                 , "Restrictions on CSOs, V-dem (sum of net changes, 5yrs)"
                                 , "Vertical accountability (mean, 5yrs)"
                                 , "Horizontal accountability (mean, 5yrs)"
                                 , "Conflict years (mean, 5yrs)"
                                 , "GDP (mean, 5yrs, logged)"
                                 , "Population (mean, 5yrs)"
                                 , "Populist leader years (count, 5yrs)"
                                 , "Any populist leader (5yrs)") )



#### Multivariate models in main body of the paper ####

#### Models 1: Public Goods Spending ####
set.seed(1234)

# DV: VDEM: Public goods spending (higher values = more public goods)

formula1a <- as.formula( paste0("VDEM_publicgoodsspend ~ RESTRICT_COUNT_NOBAN_netl5", controls,  " | cown  + year | 0 | cown") )

# Model SG1
model1a <- felm(formula = formula1a,
                data = subset(data, year>=2000 & year<=2016 & BANNED==0)  )
summary(model1a)

# Formula
formula1b <- as.formula( paste0("VDEM_publicgoodsspend ~ VDEM_CSOrepress_netl5", controls,  " | cown  + year | 0 | cown + year") )

# Model SG1
model1b <- felm(formula = formula1b,
                data = data)
summary(model1b)

stargazer(model1a, model1b
          ,covariate.labels = c("Restrictions on CSOs (sum of net changes, 5yrs)"
                                , "Restrictions on CSOs, V-dem (sum of net changes, 5yrs)"
                                , "Conflict years (mean, 5yrs)"
                                , "Vertical accountability (mean, 5yrs)"
                                , "Horizontal accountability (mean, 5yrs)"
                                , "GDP (mean, 5yrs)"
                                , "Population (mean, 5yrs)")
          , type="text"
          , model.numbers=F
          , out="Table2.txt"
          , title="Two-fixed effects linear models of public goods orientation in national spending"
          , column.labels=c("Model 1a", "Model 1b")
          , dep.var.labels=c("Public goods orientation in spending (V-Dem)")
          , notes = "(cluster-robust standard errors)")

# Interpretation: 
d <- subset(data, year>=2000 & year<=2016 & BANNED==0) 

paste0("One standard deviation increase in restrictions leads to ", 
       model1a$coefficients[1]*(0 + sd(na.omit(d$RESTRICT_COUNT_NOBAN_netl5))),
       " units increase in public goods spending, which corresponds to ",
       model1a$coefficients[1]*(0 + sd(na.omit(d$RESTRICT_COUNT_NOBAN_netl5))) / sd(na.omit(d$VDEM_publicgoodsspend)) * 100,
       " percent of its standard deviation.")

paste0("One standard deviation increase in restrictions leads to ", 
       model1b$coefficients[1]*(0 + sd(na.omit(data$VDEM_CSOrepress_netl5))),
       " units increase in public goods spending, which corresponds to ",
       model1b$coefficients[1]*(0 + sd(na.omit(data$VDEM_CSOrepress_netl5))) / sd(na.omit(data$VDEM_publicgoodsspend)) * 100,
       " percent of its standard deviation.")

paste0("One standard deviation increase in restrictions leads to ", 
       model1b$coefficients[4]*(0 + sd(na.omit(data$Control_VDEM_HorizontalAccount_netl5))),
       " units increase in public goods spending, which corresponds to ",
       model1b$coefficients[4]*(0 + sd(na.omit(data$Control_VDEM_HorizontalAccount_netl5))) / sd(na.omit(data$VDEM_publicgoodsspend)) * 100,
       " percent of its standard deviation.")



#### Models 2: Corruption ####
# DV: VDEM: Corruption (higher values = more corruption)

formula2a <- as.formula( paste0("VDEM_corruptionex ~ RESTRICT_COUNT_NOBAN_netl5", controls,  " | cown  + year | 0 | cown") )

# Model SG1
model2a <- felm(formula = formula2a,
                data = subset(data, year>=2000 & year<=2016 & BANNED==0) )
summary(model2a)


# Formula
formula2b <- as.formula( paste0("VDEM_corruptionex ~ VDEM_CSOrepress_netl5", controls,  " | cown  + year | 0 | cown") )

# Model SG1
model2b <- felm(formula = formula2b,
                data = data)
summary(model2b)


stargazer(model2a, model2b
          ,covariate.labels = c("Restrictions on CSOs (sum of net changes, 5yrs)"
                                , "Restrictions on CSOs, V-dem (sum of net changes, 5yrs)"
                                , "Conflict years (mean, 5yrs)"
                                , "Vertical accountability (mean, 5yrs)"
                                , "Horizontal accountability (mean, 5yrs)"
                                , "GDP (mean, 5yrs)"
                                , "Population (mean, 5yrs)")
          , type="text"
          , out="Table3.txt"
          , title="Two-fixed effects linear models of scale of corruption"
          , column.labels=c("Model 2a", "Model 2b")
          , dep.var.labels=c("Executive corruption index (V-Dem)")
          , notes = "(cluster-robust standard errors)")


# Interpretation: 
d <- subset(data, year>=2000 & year<=2016 & BANNED==0) 

paste0("One standard deviation increase in restrictions leads to ", 
       model2a$coefficients[1]*(0 + sd(na.omit(d$RESTRICT_COUNT_NOBAN_netl5))),
       " units increase in executive corruption, which corresponds to ",
       model2a$coefficients[1]*(0 + sd(na.omit(d$RESTRICT_COUNT_NOBAN_netl5))) / sd(na.omit(d$VDEM_corruptionex)) * 100,
       " percent of its standard deviation.")

paste0("One standard deviation increase in restrictions leads to ", 
       model2b$coefficients[1]*(0 + sd(na.omit(data$VDEM_CSOrepress_netl5))),
       " units increase in executive corruption, which corresponds to ",
       model2b$coefficients[1]*(0 + sd(na.omit(data$VDEM_CSOrepress_netl5))) / sd(na.omit(data$VDEM_corruptionex)) * 100,
       " percent of its standard deviation.")

paste0("One standard deviation increase in restrictions leads to ", 
       model2b$coefficients[4]*(0 + sd(na.omit(data$Control_VDEM_HorizontalAccount_netl5))),
       " units increase in public goods spending, which corresponds to ",
       model2b$coefficients[4]*(0 + sd(na.omit(data$Control_VDEM_HorizontalAccount_netl5))) / sd(na.omit(data$VDEM_corruptionex)) * 100,
       " percent of its standard deviation.")



#### Models 3: Clientelism: Vote buying and clientelist parties####
# DV: VDEM: Vote buying (higher values = more vote buying; scale reversed compared to V-Dem Original)

formula3a <- as.formula( paste0("VDEM_votebuy ~ RESTRICT_COUNT_NOBAN_netl5", controls,  " | cown  + year | 0 | cown") )

# Model SG1
model3a <- felm(formula = formula3a,
                data = subset(data, year>=2000 & year<=2016 & BANNED==0) )
summary(model3a)

# Formula
formula3b <- as.formula( paste0("VDEM_votebuy ~ VDEM_CSOrepress_netl5", controls,  " | cown  + year | 0 | cown") )

# Model SG1
model3b <- felm(formula = formula3b,
                data = data)
summary(model3b)


# Models 3: Clientelism: Clientelist parties
# DV: VDEM: Clientelist parties (higher values = more clientelist parties; scale reversed compared to V-Dem Original)

formula3c <- as.formula( paste0("VDEM_progrClientParties ~ RESTRICT_COUNT_NOBAN_netl5", controls,  " | cown  + year | 0 | cown") )

# Model SG1
model3c <- felm(formula = formula3c,
                data = subset(data, year>=2000 & year<=2016 & BANNED==0) )
summary(model3c)

# Formula
formula3d <- as.formula( paste0("VDEM_progrClientParties ~ VDEM_CSOrepress_netl5", controls,  " | cown  + year | 0 | cown") )

# Model SG1
model3d <- felm(formula = formula3d,
                data = data)
summary(model3d)


stargazer(model3a, model3b, model3c, model3d
          ,covariate.labels = c("Restrictions on CSOs (sum of net changes, 5yrs)"
                                , "Restrictions on CSOs, V-dem (sum of net changes, 5yrs)"
                                , "Conflict years (mean, 5yrs)"
                                , "Vertical accountability (mean, 5yrs)"
                                , "Horizontal accountability (mean, 5yrs)"
                                , "GDP (mean, 5yrs)"
                                , "Population (mean, 5yrs)")
          , type="text"
          , out="Table4.txt"
          , title="Two-fixed effects linear models of scale of clientelism"
          , column.labels=c("Model 3a", "Model 3b","Model 3c", "Model 3d")
          , dep.var.labels=c("Vote buying","Clientelist Parties")
          , notes = "(cluster-robust standard errors)")


# Interpretation: 
d <- subset(data, year>=2000 & year<=2016 & BANNED==0) 

paste0("One standard deviation increase in restrictions leads to ", 
       model3a$coefficients[1]*(0 + sd(na.omit(d$RESTRICT_COUNT_NOBAN_netl5))),
       " units increase in vote buying, which corresponds to ",
       model3a$coefficients[1]*(0 + sd(na.omit(d$RESTRICT_COUNT_NOBAN_netl5))) / sd(na.omit(d$VDEM_votebuy)) * 100,
       " percent of its standard deviation.")

paste0("One standard deviation increase in restrictions leads to ", 
       model3b$coefficients[1]*(0 + sd(na.omit(data$VDEM_CSOrepress_netl5))),
       " units increase in public goods spending, which corresponds to ",
       model3b$coefficients[1]*(0 + sd(na.omit(data$VDEM_CSOrepress_netl5))) / sd(na.omit(data$VDEM_votebuy)) * 100,
       " percent of its standard deviation.")


paste0("One standard deviation increase in restrictions leads to ", 
       model3c$coefficients[1]*(0 + sd(na.omit(d$RESTRICT_COUNT_NOBAN_netl5))),
       " units increase in clientelism, which corresponds to ",
       model3c$coefficients[1]*(0 + sd(na.omit(d$RESTRICT_COUNT_NOBAN_netl5))) / sd(na.omit(d$VDEM_progrClientParties)) * 100,
       " percent of its standard deviation.")

paste0("One standard deviation increase in restrictions leads to ", 
       model3d$coefficients[1]*(0 + sd(na.omit(data$VDEM_CSOrepress_netl5))),
       " units increase in clientelism, which corresponds to ",
       model3d$coefficients[1]*(0 + sd(na.omit(data$VDEM_CSOrepress_netl5))) / sd(na.omit(data$VDEM_progrClientParties)) * 100,
       " percent of its standard deviation.")


paste0("One standard deviation increase in horizontal accountability leads to ", 
       model3a$coefficients[4]*(0 + sd(na.omit(d$Control_VDEM_HorizontalAccount_netl5))),
       " units increase in vote buying, which corresponds to ",
       model3a$coefficients[4]*(0 + sd(na.omit(d$Control_VDEM_HorizontalAccount_netl5))) / sd(na.omit(d$VDEM_votebuy)) * 100,
       " percent of its standard deviation.")

paste0("One standard deviation increase in horizontal accountability leads to ", 
       model3b$coefficients[4]*(0 + sd(na.omit(data$Control_VDEM_HorizontalAccount_netl5))),
       " units increase in vote buying, which corresponds to ",
       model3b$coefficients[4]*(0 + sd(na.omit(data$Control_VDEM_HorizontalAccount_netl5))) / sd(na.omit(data$VDEM_votebuy)) * 100,
       " percent of its standard deviation.")

paste0("One standard deviation increase in horizontal accountabilityleads to ", 
       model3c$coefficients[4]*(0 + sd(na.omit(d$Control_VDEM_HorizontalAccount_netl5))),
       " units increase in clientelism, which corresponds to ",
       model3c$coefficients[4]*(0 + sd(na.omit(d$Control_VDEM_HorizontalAccount_netl5))) / sd(na.omit(d$VDEM_progrClientParties)) * 100,
       " percent of its standard deviation.")

paste0("One standard deviation increase in horizontal accountability leads to ", 
       model3d$coefficients[4]*(0 + sd(na.omit(data$Control_VDEM_HorizontalAccount_netl5))),
       " units increase in clientelism, which corresponds to ",
       model3d$coefficients[4]*(0 + sd(na.omit(data$Control_VDEM_HorizontalAccount_netl5))) / sd(na.omit(data$VDEM_progrClientParties)) * 100,
       " percent of its standard deviation.")







#### Test of mechanisms #####


#### Models 4: Engaged society ####
# DV: VDEM: (higher values = more engagement)

formula4a <- as.formula( paste0("Mech_VDEM_engagedSociety_netl5 ~ RESTRICT_COUNT_NOBAN_netl5", controls,  " | cown  + year | 0 | cown") )

# Model SG1
model4a <- felm(formula = formula4a,
                data = subset(data, year>=2000 & year<=2016 & BANNED==0) )
summary(model4a)

# Substantive interpretation
d <- subset(data, year>=2000 & year<=2016 & BANNED==0)
(model4a$coefficients[1]*sd(na.omit(d$RESTRICT_COUNT_NOBAN_netl5)))/ sd(d$Mech_VDEM_engagedSociety_netl5) * 100

# Formula
formula4b <- as.formula( paste0("Mech_VDEM_engagedSociety_netl5 ~ VDEM_CSOrepress_netl5", controls,  " | cown  + year | 0 | cown") )

# Model SG1
model4b <- felm(formula = formula4b,
                data = data)
summary(model4b)

# Substantive interpretation
model4b$coefficients[1]*sd(na.omit(data$VDEM_CSOrepress_netl5)) / sd(data$Mech_VDEM_engagedSociety_netl5) * 100



stargazer(model4a, model4b
          ,covariate.labels = c("Restrictions on CSOs (sum of net changes, 5yrs)"
                                , "Restrictions on CSOs, V-dem (sum of net changes, 5yrs)"
                                , "Conflict years (mean, 5yrs)"
                                , "Conflict years (mean, 5yrs)"
                                , "Vertical accountability (mean, 5yrs)"
                                , "Horizontal accountability (mean, 5yrs)"
                                , "GDP (mean, 5yrs)"
                                , "Population (mean, 5yrs)")
          , type="text"
          , out="Table5.txt"
          , title="Two-fixed effects linear models of engaged society"
          , column.labels=c("Model 4a", "Model 4b")
          , dep.var.labels=c("Engaged society")
          , notes = "(cluster-robust standard errors)")


#### Model 5: Naming & Shaming ####
# DV: HRO conflictual events targeting government (Murdie & Davis)
data$shamingINGO
formula5a <- as.formula( paste0("Mech_shamingINGO_netl5 ~ RESTRICT_COUNT_NOBAN_netl5", controls,  " | cown  + year | 0 | cown") )

# Model SG1
model5a <- felm(formula = formula5a,
                data = subset(data, year>=2000 & year<=2007 & BANNED==0) )
summary(model5a)


# Formula
formula5b <- as.formula( paste0("Mech_shamingINGO_netl5 ~ VDEM_CSOrepress_netl5", controls,  " | cown  + year | 0 | cown") )

# Model SG1
model5b <- felm(formula = formula5b,
                data = subset(data, year<=2007) )
summary(model5b)


stargazer(model5a, model5b
          ,covariate.labels = c("Restrictions on CSOs (sum of net changes, 5yrs)"
                                , "Restrictions on CSOs, V-dem (sum of net changes, 5yrs)"
                                , "Conflict years (mean, 5yrs)"
                                , "Conflict years (mean, 5yrs)"
                                , "Vertical accountability (mean, 5yrs)"
                                , "Horizontal accountability (mean, 5yrs)"
                                , "GDP (mean, 5yrs)"
                                , "Population (mean, 5yrs)")
          , type="text"
          , out="Table6.txt"
          , title="Two-fixed effects linear models of international shaming"
          , column.labels=c("Model 5a", "Model 5b")
          , dep.var.labels=c("Shaming (Murdie and Davis)")
          , notes = "(cluster-robust standard errors)")



#### Models 6: Protest ####
# DV: Clark and Regan Protest Count (binary)

formula6a <- as.formula( paste0("Mech_CR_protest_netl5_log ~ RESTRICT_COUNT_NOBAN_netl5", controls,  " | cown  + year | 0 | cown") )

# Model SG1
model6a <- felm(formula = formula6a,
                data = subset(data, year>=2000 & year<=2016 & BANNED==0) )
summary(model6a)


# Formula
formula6b <- as.formula( paste0("Mech_CR_protest_netl5_log ~ VDEM_CSOrepress_netl5", controls,  " | cown  + year | 0 | cown") )

# Model SG1
model6b <- felm(formula = formula6b,
                data = subset(data, year<=2016) )
summary(model6b)

#Interpretation
(model6b$coefficients[1]* sd(na.omit(data$VDEM_CSOrepress_netl5))) /sd(na.omit(data$Mech_CR_protest_netl5_log))*100


stargazer(model6a, model6b
          ,covariate.labels = c("Restrictions on CSOs (sum of net changes, 5yrs)"
                                , "Restrictions on CSOs, V-dem (sum of net changes, 5yrs)"
                                , "Conflict years (mean, 5yrs)"
                                , "Conflict years (mean, 5yrs)"
                                , "Vertical accountability (mean, 5yrs)"
                                , "Horizontal accountability (mean, 5yrs)"
                                , "GDP (mean, 5yrs)"
                                , "Population (mean, 5yrs)")
          , type="text"
          , out="Table7.txt"
          , title="Two-fixed effects linear models of protest"
          , column.labels=c("Model 6a", "Model 6b")
          , dep.var.labels=c("Protest (Clark and Regan)")
          , notes = "(cluster-robust standard errors)")







#### Predictions (out of sample) #####

controls <- " + Control_UCDP_ongoing_conflict_netl5 +  Control_WB_controls_GDP_log_netl5 + Control_WB_controls_population_log_netl5"

## Predictive power of restrictions from Bakke et al.
formula1a_base <- as.formula( paste0("VDEM_publicgoodsspend  ~ ", controls, "+as.factor(cown) +as.factor(year)") )
formula1a_restrict <- as.formula( paste0("VDEM_publicgoodsspend  ~  RESTRICT_COUNT_NOBAN_netl5", controls, "+as.factor(cown) +as.factor(year)") )
formula1a_vertic <- as.formula( paste0("VDEM_publicgoodsspend ~ Control_VDEM_VerticalAccount_netl5_pred", controls, "+as.factor(cown) +as.factor(year)") )
formula1a_horiz <- as.formula( paste0("VDEM_publicgoodsspend ~ Control_VDEM_HorizontalAccount_netl5_pred", controls, "+as.factor(cown) +as.factor(year)") )

myControl1 = trainControl(method = "cv", number=4)
myControl2 = trainControl(method="repeatedcv", number=4, repeats=5)

model_base <- train(formula1a_base,
                    data = subset(data, year>=2000 & year<=2016 & BANNED==0),                     
                    trControl = myControl2,             
                    method = "lm",                      
                    na.action = na.pass)  

model_restrict <- train(formula1a_restrict,
                        data = subset(data, year>=2000 & year<=2016 & BANNED==0),                      
                        trControl = myControl2,             
                        method = "lm",                      
                        na.action = na.pass)  

model_vertic <- train(formula1a_vertic,
                      data = subset(data, year>=2000 & year<=2016 & BANNED==0),                      
                      trControl = myControl2,             
                      method = "lm",                      
                      na.action = na.pass)  

model_horiz <- train(formula1a_horiz,
                     data = subset(data, year>=2000 & year<=2016 & BANNED==0),                      
                     trControl = myControl2,             
                     method = "lm",                      
                     na.action = na.pass)  

# Results
model_base
model_restrict
model_vertic
model_horiz

result_base<-model_base$results$Rsquared
result_restrict<-model_restrict$results$Rsquared
result_vertic<-model_vertic$results$Rsquared
result_horiz<-model_horiz$results$Rsquared

results1 <- c(result_restrict,result_vertic,result_horiz)-result_base
results1 <- as.data.frame(cbind(results1, c("Net changes in restrictions\nSample: Bakke et al. (1996-2016)", "Net changes in vertical account.\nSample: Bakke et al. (1996-2016)"
                                            , "Net changes in horizontal account.\nSample: Bakke et al. (1996-2016)")))

## Predictive power of V-Dem restrictions
formula1b_base <- as.formula( paste0("VDEM_publicgoodsspend ~ ", controls, "+as.factor(cown) +as.factor(year)") )
formula1b_restrict <- as.formula( paste0("VDEM_publicgoodsspend ~ VDEM_CSOrepress_netl5", controls, "+as.factor(cown) +as.factor(year)") )
formula1b_vertic <- as.formula( paste0("VDEM_publicgoodsspend ~ Control_VDEM_VerticalAccount_netl5_pred", controls, "+as.factor(cown) +as.factor(year)") )
formula1b_horiz <- as.formula( paste0("VDEM_publicgoodsspend ~ Control_VDEM_HorizontalAccount_netl5_pred", controls, "+as.factor(cown) +as.factor(year)") )

model_base <- train(formula1b_base,
                    data = data,                        
                    trControl = myControl2,             
                    method = "lm",                      
                    na.action = na.pass)  

model_restrict <- train(formula1b_restrict,
                        data = data,                        
                        trControl = myControl2,             
                        method = "lm",                      
                        na.action = na.pass)  

model_vertic <- train(formula1b_vertic,
                      data = data,                        
                      trControl = myControl2,             
                      method = "lm",                      
                      na.action = na.pass)  

model_horiz <- train(formula1b_horiz,
                     data = data,                        
                     trControl = myControl2,             
                     method = "lm",                      
                     na.action = na.pass)  
# Results
model_base
model_restrict
model_vertic
model_horiz


result_base<-model_base$results$Rsquared
result_restrict<-model_restrict$results$Rsquared
result_vertic<-model_vertic$results$Rsquared
result_horiz<-model_horiz$results$Rsquared
results2 <- c(result_restrict,result_vertic,result_horiz)-result_base
results2 <- as.data.frame(cbind(results2, c("Net change in restrictions\nSample: V-dem (1990-2021)", "Net change in vertical account.\nSample: V-dem (1990-2021)"
                                            , "Net change in horizontal account.\nSample: V-dem (1990-2021)")))

# Combine results for both Bakke et al and V-Dem
colnames(results1) <- c("Rsq", "Name")
colnames(results2) <- c("Rsq", "Name")
results <- rbind(results1, results2)
results$Rsq<-as.numeric(results$Rsq)*100
results$restrictions <- rep(c("Restrictions","Institutions","Institutions"),2)

#### Figure 6: Predictive power of restrictions ... ####

png("Fig6.png", res=100)
ggbarplot(results[c(2,3,5,6,1,4),], x = "Name", y = "Rsq",
          fill = "restrictions",           # change fill color by mpg_level
          color = "white",            # Set bar border colors to white
          palette = c("grey60", "#FF3333"),            # jco journal color palett. see ?ggpar
          sort.by.groups = FALSE,     # Don't sort inside each group
          x.txtt.angle = 90,          # Rotate vertically x axis texts
          ylab = "Percentage-points of R-squared added \n compared to base model, cross-validated (k=4)",
          xlab="",
          legend.title = "",
          rotate = T,
          cex.txtt=1,
          ggtheme = theme_minimal() )
dev.off()






#### APPENDIX ####

### Appendix A ####

#### Table A1: VDEM: Public goods spending (higher values = more public goods) ####

formulaA1a <- as.formula( paste0("VDEM_publicgoodsspend ~ RESTRICT_COUNT_NOBAN_netl2", controls,  " | cown  + year | 0 | cown") )
formulaA1b <- as.formula( paste0("VDEM_publicgoodsspend ~ RESTRICT_COUNT_NOBAN_netl8", controls,  " | cown  + year | 0 | cown") )

# Model SG1
modelA1a <- felm(formula = formulaA1a,
                 data = subset(data, year>=2000 & year<=2016 & BANNED==0)  )
summary(modelA1a)

modelA1b <- felm(formula = formulaA1b,
                 data = subset(data, year>=2000 & year<=2016 & BANNED==0)  )
summary(modelA1b)

# Formula
formulaA1c <- as.formula( paste0("VDEM_publicgoodsspend ~ VDEM_CSOrepress_netl2", controls,  " | cown  + year | 0 | cown + year") )
formulaA1d <- as.formula( paste0("VDEM_publicgoodsspend ~ VDEM_CSOrepress_netl8", controls,  " | cown  + year | 0 | cown + year") )

# Model SG1
modelA1c <- felm(formula = formulaA1c,
                 data = data)
summary(modelA1c)

modelA1d <- felm(formula = formulaA1d,
                 data = data)
summary(modelA1d)

stargazer(modelA1a, modelA1b, modelA1c, modelA1d
          ,covariate.labels = c("Restrictions on CSOs (sum of net changes, 2yrs)"
                                , "Restrictions on CSOs (sum of net changes, 8yrs)"
                                , "Restrictions on CSOs, V-dem (sum of net changes, 2yrs)"
                                , "Restrictions on CSOs, V-dem (sum of net changes, 8yrs)"
                                , "Conflict years (mean, 5yrs)"
                                , "Vertical accountability (mean, 5yrs)"
                                , "Horizontal accountability (mean, 5yrs)"
                                , "GDP (mean, 5yrs)"
                                , "Population (mean, 5yrs)")
          , type="text"
          , model.numbers=F
          , out="App_TableA1_PublicGoodsSpend.txt"
          , title="Two-fixed effects linear models of public goods spending"
          , column.labels=c("Model A1a", "Model A1b", "Model A1c", "Model A1d")
          , dep.var.labels=c("Public goods orientation in national spending (V-Dem)")
          , notes = "(cluster-robust standard errors)")


#### Table A2: VDEM: Corruption (higher values = more corruption) ####

formulaA2a <- as.formula( paste0("VDEM_corruptionex ~ RESTRICT_COUNT_NOBAN_netl2", controls,  " | cown  + year | 0 | cown") )
formulaA2b <- as.formula( paste0("VDEM_corruptionex ~ RESTRICT_COUNT_NOBAN_netl8", controls,  " | cown  + year | 0 | cown") )

# Model SG1
modelA2a <- felm(formula = formulaA2a,
                 data = subset(data, year>=2000 & year<=2016 & BANNED==0) )
summary(modelA2a)

modelA2b <- felm(formula = formulaA2b,
                 data = subset(data, year>=2000 & year<=2016 & BANNED==0) )
summary(modelA2b)


# Formula
formulaA2c <- as.formula( paste0("VDEM_corruptionex ~ VDEM_CSOrepress_netl2", controls,  " | cown  + year | 0 | cown") )
formulaA2d <- as.formula( paste0("VDEM_corruptionex ~ VDEM_CSOrepress_netl8", controls,  " | cown  + year | 0 | cown") )

# Model SG1
modelA2c <- felm(formula = formulaA2c,
                 data = data)
summary(modelA2c)

modelA2d <- felm(formula = formulaA2d,
                 data = data)
summary(modelA2d)


stargazer(modelA2a, modelA2b, modelA2c, modelA2d
          ,covariate.labels = c("Restrictions on CSOs (sum of net changes, 2yrs)"
                                ,"Restrictions on CSOs (sum of net changes, 8yrs)"
                                , "Restrictions on CSOs, V-dem (sum of net changes, 2yrs)"
                                , "Restrictions on CSOs, V-dem (sum of net changes, 8yrs)"
                                , "Conflict years (mean, 5yrs)"
                                , "Vertical accountability (mean, 5yrs)"
                                , "Horizontal accountability (mean, 5yrs)"
                                , "GDP (mean, 5yrs)"
                                , "Population (mean, 5yrs)")
          , type="text"
          , model.numbers=F
          , out="App_TableA2_CorruptionIndex.txt"
          , title="Two-fixed effects linear models of scale of corruption"
          , column.labels=c("Model A2a", "Model A2b", "Model A2c", "Model A2d")
          , dep.var.labels=c("Executive corruption index (V-Dem)")
          , notes = "(cluster-robust standard errors)")



#### Table A3: VDEM: Vote buying ####

# Models 3: Clientelism: Vote buying
# DV: VDEM: Vote buying (higher values = more vote buying; scale reversed compared to V-Dem Original)

formulaA3a <- as.formula( paste0("VDEM_votebuy ~ RESTRICT_COUNT_NOBAN_netl2", controls,  " | cown  + year | 0 | cown") )
formulaA3b <- as.formula( paste0("VDEM_votebuy ~ RESTRICT_COUNT_NOBAN_netl8", controls,  " | cown  + year | 0 | cown") )

# Model SG1
modelA3a <- felm(formula = formulaA3a,
                 data = subset(data, year>=2000 & year<=2016 & BANNED==0) )
summary(modelA3a)

modelA3b <- felm(formula = formulaA3b,
                 data = subset(data, year>=2000 & year<=2016 & BANNED==0) )
summary(modelA3b)

# Formula
formulaA3c <- as.formula( paste0("VDEM_votebuy ~ VDEM_CSOrepress_netl2", controls,  " | cown  + year | 0 | cown") )
formulaA3d <- as.formula( paste0("VDEM_votebuy ~ VDEM_CSOrepress_netl8", controls,  " | cown  + year | 0 | cown") )

# Model SG1
modelA3c <- felm(formula = formulaA3c,
                 data = data)
summary(modelA3c)

modelA3d <- felm(formula = formulaA3d,
                 data = data)
summary(modelA3d)


stargazer(modelA3a, modelA3b, modelA3c, modelA3d
          ,covariate.labels = c("Restrictions on CSOs (sum of net changes, 2yrs)"
                                ,"Restrictions on CSOs (sum of net changes, 8yrs)"
                                , "Restrictions on CSOs, V-dem (sum of net changes, 2yrs)"
                                , "Restrictions on CSOs, V-dem (sum of net changes, 8yrs)"
                                , "Conflict years (mean, 5yrs)"
                                , "Vertical accountability (mean, 5yrs)"
                                , "Horizontal accountability (mean, 5yrs)"
                                , "GDP (mean, 5yrs)"
                                , "Population (mean, 5yrs)")
          , type="text"
          , model.numbers=F
          , out="App_TableA3_Clientelism_Vote buying.txt"
          , title="Two-fixed effects linear models of scale of clientelism"
          , column.labels=c("Model A3a", "Model A3b","Model A3c", "Model A3d")
          , dep.var.labels=c("Vote buying")
          , notes = "(cluster-robust standard errors)")


#### Table A4: VDEM: Clientelist parties ####

# Models 4: Clientelism: Clientelist parties
# DV: VDEM: Clientelist parties (higher values = more clientelist parties; scale reversed compared to V-Dem Original)

formulaA4a <- as.formula( paste0("VDEM_progrClientParties ~ RESTRICT_COUNT_NOBAN_netl2", controls,  " | cown  + year | 0 | cown") )
formulaA4b <- as.formula( paste0("VDEM_progrClientParties ~ RESTRICT_COUNT_NOBAN_netl8", controls,  " | cown  + year | 0 | cown") )

# Model SG1
modelA4a <- felm(formula = formulaA4a,
                 data = subset(data, year>=2000 & year<=2016 & BANNED==0) )
summary(modelA4a)

modelA4b <- felm(formula = formulaA4b,
                 data = subset(data, year>=2000 & year<=2016 & BANNED==0) )
summary(modelA4b)

# Formula
formulaA4c <- as.formula( paste0("VDEM_progrClientParties ~ VDEM_CSOrepress_netl2", controls,  " | cown  + year | 0 | cown") )
formulaA4d <- as.formula( paste0("VDEM_progrClientParties ~ VDEM_CSOrepress_netl8", controls,  " | cown  + year | 0 | cown") )

# Model SG1
modelA4c <- felm(formula = formulaA4c,
                 data = data)
summary(modelA4c)

modelA4d <- felm(formula = formulaA4d,
                 data = data)
summary(modelA4d)



stargazer( modelA4a, modelA4b, modelA4c, modelA4d
           ,covariate.labels = c("Restrictions on CSOs (sum of net changes, 2yrs)"
                                 ,"Restrictions on CSOs (sum of net changes, 8yrs)"
                                 , "Restrictions on CSOs, V-dem (sum of net changes, 2yrs)"
                                 , "Restrictions on CSOs, V-dem (sum of net changes, 8yrs)"
                                 , "Conflict years (mean, 5yrs)"
                                 , "Vertical accountability (mean, 5yrs)"
                                 , "Horizontal accountability (mean, 5yrs)"
                                 , "GDP (mean, 5yrs)"
                                 , "Population (mean, 5yrs)")
           , type="text"
           , model.numbers=F
           , out="App_TableA4_Clientelism_Parties.txt"
           , title="Two-fixed effects linear models of scale of clientelism"
           , column.labels=c("Model A4a", "Model A4b","Model A4c", "Model A4d")
           , dep.var.labels=c("Clientelist Parties")
           , notes = "(cluster-robust standard errors)")




#### Appendix B ####

#### Table B1 ####

# Control variables
controls2 <- paste0(" + Control_UCDP_ongoing_conflict_netl5 + Control_VDEM_VerticalAccount_netl5 + Control_VDEM_HorizontalAccount_netl5 + Control_WB_controls_GDP_log_netl5 + Control_WB_controls_population_log_netl5 + Control_populistleader_l5yrs_bi")

# DV: VDEM: Public goods spending (higher values = more public goods)

formulaB1b <- as.formula( paste0("VDEM_publicgoodsspend ~ RESTRICT_COUNT_NOBAN_netl5", controls2,  " | cown  + year | 0 | cown") )

# Model SG1


modelB1b <- felm(formula = formulaB1b,
                 data = subset(data, year>=2000 & year<=2016 & BANNED==0)  )
summary(modelB1b)

# Formula
formulaB1d <- as.formula( paste0("VDEM_publicgoodsspend ~ VDEM_CSOrepress_netl5", controls2,  " | cown  + year | 0 | cown + year") )

# Model SG1

modelB1d <- felm(formula = formulaB1d,
                 data = data)
summary(modelB1d)

stargazer(modelB1b, modelB1d
          ,covariate.labels = c("Restrictions on CSOs (sum of net changes, 5yrs)"
                                , "Restrictions on CSOs, V-dem (sum of net changes, 5yrs)"
                                , "Conflict years (mean, 5yrs)"
                                , "Vertical accountability (mean, 5yrs)"
                                , "Horizontal accountability (mean, 5yrs)"
                                , "GDP (mean, 5yrs)"
                                , "Population (mean, 5yrs)"
                                , "Any populist leader (5yrs)"
                                , "Any populist leader (5yrs)")
          , type="text"
          , model.numbers=F
          , out="App_TableB1_PublicGoodsSpend.txt"
          , title="Two-fixed effects linear models of public goods spending (with populist leader spells as control)"
          , column.labels=c("Model B1a", "Model B1b")
          , dep.var.labels=c("Public goods orientation in national spending (V-Dem)")
          , notes = "(cluster-robust standard errors)")


#### Table B2 ####

# Models 2: Corruption
# DV: VDEM: Corruption (higher values = more corruption)

formulaB2b <- as.formula( paste0("VDEM_corruptionex ~ RESTRICT_COUNT_NOBAN_netl5", controls2,  " | cown  + year | 0 | cown") )

# Model SG1

modelB2b <- felm(formula = formulaB2b,
                 data = subset(data, year>=2000 & year<=2016 & BANNED==0) )
summary(modelB2b)


# Formula
formulaB2d <- as.formula( paste0("VDEM_corruptionex ~ VDEM_CSOrepress_netl5", controls2,  " | cown  + year | 0 | cown") )

# Model SG1
modelB2d <- felm(formula = formulaB2d,
                 data = data)
summary(modelB2d)


stargazer(modelB2b, modelB2d
          ,covariate.labels = c("Restrictions on CSOs (sum of net changes, 5yrs)"
                                , "Restrictions on CSOs, V-dem (sum of net changes, 5yrs)"
                                , "Conflict years (mean, 5yrs)"
                                , "Vertical accountability (mean, 5yrs)"
                                , "Horizontal accountability (mean, 5yrs)"
                                , "GDP (mean, 5yrs)"
                                , "Population (mean, 5yrs)"
                                , "Populist leader years (count, 5yrs)"
                                , "Any populist leader (5yrs)")
          , type="text"
          , model.numbers=F
          , out="App_TableB2_CorruptionIndex.txt"
          , title="Two-fixed effects linear models of scale of corruption  (with populist leader spells as control)"
          , column.labels=c("Model B2a", "Model B2b")
          , dep.var.labels=c("Executive corruption index (V-Dem)")
          , notes = "(cluster-robust standard errors)")



#### Appendix C: Level of Restrictions on CSOs ####

#### Table C1 ####

set.seed(1234)

# DV: VDEM: Public goods spending (higher values = more public goods)
formula1a <- as.formula( paste0("VDEM_publicgoodsspend ~ RESTRICT_COUNT_NOBAN_netl5 + RESTRICT_COUNT_NOBAN_nochange_netl5", controls,  " | cown  + year | 0 | cown") )
formula1b <- as.formula( paste0("VDEM_publicgoodsspend ~ RESTRICT_COUNT_NOBAN_nochange_netl5", controls,  " | cown  + year | 0 | cown") )

# Model SG1
model1a <- felm(formula = formula1a,
                data = subset(data, year>=2000 & year<=2016 & BANNED==0)  )
summary(model1a)

model1b <- felm(formula = formula1b,
                data = subset(data, year>=2000 & year<=2016 & BANNED==0)  )
summary(model1b)

# Formula
formula1c <- as.formula( paste0("VDEM_publicgoodsspend ~ VDEM_CSOrepress_netl5 + VDEM_CSOrepress_nochange_netl5", controls,  " | cown  + year | 0 | cown + year") )
formula1d <- as.formula( paste0("VDEM_publicgoodsspend ~ VDEM_CSOrepress_nochange_netl5", controls,  " | cown  + year | 0 | cown + year") )

# Model SG1
model1c <- felm(formula = formula1c,
                data = data)
summary(model1c)

model1d <- felm(formula = formula1d,
                data = data)
summary(model1d)


stargazer(model1a, model1b,model1c, model1d
          ,covariate.labels = c("Restrictions on CSOs (sum of net changes, 5yrs)"
                                ,"Restrictions on CSOs (level, 5yrs)"
                                , "Restrictions on CSOs, V-dem (sum of net changes, 5yrs)"
                                , "Restrictions on CSOs, V-dem (level, 5yrs)"
                                , "Vertical accountability (mean, 5yrs)"
                                , "Horizontal accountability (mean, 5yrs)"
                                , "Conflict years (mean, 5yrs)" 
                                , "GDP (mean, 5yrs)"
                                , "Population (mean, 5yrs)")
          , type="text"
          , model.numbers=F
          , out="App_TableC1_PublicGoodsSpend_LevelRestrict.txt"
          , title="Two-fixed effects linear models of public goods spending"
          , column.labels=c("Model C1a", "Model C1b", "Model C1c", "Model C1d")
          , dep.var.labels=c("Public goods orientation in national spending (V-Dem)")
          , notes = "(cluster-robust standard errors)")


#### Table C2 ####

# Redefine Control Variables
controls_new <- paste0(" + Control_UCDP_ongoing_conflict_netl5 + Control_WB_controls_GDP_log_netl5 + Control_WB_controls_population_log_netl5")

set.seed(1234)

# DV: VDEM: Public goods spending (higher values = more public goods)
formula1a <- as.formula( paste0("VDEM_publicgoodsspend ~ RESTRICT_COUNT_NOBAN_netl5 + RESTRICT_COUNT_NOBAN_nochange_netl5", controls_new,  " | cown  + year | 0 | cown") )
formula1b <- as.formula( paste0("VDEM_publicgoodsspend ~ RESTRICT_COUNT_NOBAN_nochange_netl5", controls_new,  " | cown  + year | 0 | cown") )

# Model SG1
model1a <- felm(formula = formula1a,
                data = subset(data, year>=2000 & year<=2016 & BANNED==0)  )
summary(model1a)

model1b <- felm(formula = formula1b,
                data = subset(data, year>=2000 & year<=2016 & BANNED==0)  )
summary(model1b)



# Formula
formula1c <- as.formula( paste0("VDEM_publicgoodsspend ~ VDEM_CSOrepress_netl5 + VDEM_CSOrepress_nochange_netl5", controls_new,  " | cown  + year | 0 | cown + year") )
formula1d <- as.formula( paste0("VDEM_publicgoodsspend ~ VDEM_CSOrepress_nochange_netl5", controls_new,  " | cown  + year | 0 | cown + year") )

# Model SG1
model1c <- felm(formula = formula1c,
                data = data)
summary(model1c)

model1d <- felm(formula = formula1d,
                data = data)
summary(model1d)

stargazer(model1a, model1b,model1c, model1d
          ,covariate.labels = c("Restrictions on CSOs (sum of net changes, 5yrs)"
                                ,"Restrictions on CSOs (level, 5yrs)"
                                , "Restrictions on CSOs, V-dem (sum of net changes, 5yrs)"
                                , "Restrictions on CSOs, V-dem (level, 5yrs)"
                                , "Conflict years (mean, 5yrs)"  
                                , "GDP (mean, 5yrs)"
                                , "Population (mean, 5yrs)")
          , type="text"
          , model.numbers=F
          , out="App_TableC2_PublicGoodsSpend_LevelRestrict_NewControls.txt"
          , title="Two-fixed effects linear models of public goods spending, different controls"
          , column.labels=c("Model C2a", "Model C2b", "Model C2c", "Model C2d")
          , dep.var.labels=c("Public goods orientation in national spending (V-Dem)")
          , notes = "(cluster-robust standard errors)")



#### Table C3 ####

# Standard Control Variables

set.seed(1234)

# DV: VDEM: Corruption
formula1a <- as.formula( paste0("VDEM_corruptionex ~ RESTRICT_COUNT_NOBAN_netl5 + RESTRICT_COUNT_NOBAN_nochange_netl5", controls,  " | cown  + year | 0 | cown") )
formula1b <- as.formula( paste0("VDEM_corruptionex ~ RESTRICT_COUNT_NOBAN_nochange_netl5", controls,  " | cown  + year | 0 | cown") )

# Model SG1
model1a <- felm(formula = formula1a,
                data = subset(data, year>=2000 & year<=2016 & BANNED==0)  )
summary(model1a)

model1b <- felm(formula = formula1b,
                data = subset(data, year>=2000 & year<=2016 & BANNED==0)  )
summary(model1b)

# Formula
formula1c <- as.formula( paste0("VDEM_corruptionex ~ VDEM_CSOrepress_netl5 + VDEM_CSOrepress_nochange_netl5", controls,  " | cown  + year | 0 | cown + year") )
formula1d <- as.formula( paste0("VDEM_corruptionex ~ VDEM_CSOrepress_nochange_netl5", controls,  " | cown  + year | 0 | cown + year") )

# Model SG1
model1c <- felm(formula = formula1c,
                data = data)
summary(model1c)

model1d <- felm(formula = formula1d,
                data = data)
summary(model1d)


stargazer(model1a, model1b,model1c, model1d
          ,covariate.labels = c("Restrictions on CSOs (sum of net changes, 5yrs)"
                                ,"Restrictions on CSOs (level, 5yrs)"
                                , "Restrictions on CSOs, V-dem (sum of net changes, 5yrs)"
                                , "Restrictions on CSOs, V-dem (level, 5yrs)"
                                , "Vertical accountability (mean, 5yrs)"
                                , "Horizontal accountability (mean, 5yrs)"
                                , "Conflict years (mean, 5yrs)" 
                                , "GDP (mean, 5yrs)"
                                , "Population (mean, 5yrs)")
          , type="text"
          , model.numbers=F
          , out="App_TableC3_Corruption_LevelRestrict.txt"
          , title="Two-fixed effects linear models of executive corruption"
          , column.labels=c("Model 3a", "Model 3b", "Model 3c", "Model 3d")
          , dep.var.labels=c("Executive corruption index (V-Dem)")
          , notes = "(cluster-robust standard errors)")


#### Table C4 ####

# Redefine Control Variables
controls_new <- paste0(" + Control_UCDP_ongoing_conflict_netl5 + Control_WB_controls_GDP_log_netl5 + Control_WB_controls_population_log_netl5")

set.seed(1234)

# DV: VDEM: Corruption
formula1a <- as.formula( paste0("VDEM_corruptionex ~ RESTRICT_COUNT_NOBAN_netl5 + RESTRICT_COUNT_NOBAN_nochange_netl5", controls_new,  " | cown  + year | 0 | cown") )
formula1b <- as.formula( paste0("VDEM_corruptionex ~ RESTRICT_COUNT_NOBAN_nochange_netl5", controls_new,  " | cown  + year | 0 | cown") )

# Model SG1
model1a <- felm(formula = formula1a,
                data = subset(data, year>=2000 & year<=2016 & BANNED==0)  )
summary(model1a)

model1b <- felm(formula = formula1b,
                data = subset(data, year>=2000 & year<=2016 & BANNED==0)  )
summary(model1b)



# Formula
formula1c <- as.formula( paste0("VDEM_corruptionex ~ VDEM_CSOrepress_netl5 + VDEM_CSOrepress_nochange_netl5", controls_new,  " | cown  + year | 0 | cown + year") )
formula1d <- as.formula( paste0("VDEM_corruptionex ~ VDEM_CSOrepress_nochange_netl5", controls_new,  " | cown  + year | 0 | cown + year") )

# Model SG1
model1c <- felm(formula = formula1c,
                data = data)
summary(model1c)

model1d <- felm(formula = formula1d,
                data = data)
summary(model1d)

stargazer(model1a, model1b,model1c, model1d
          ,covariate.labels = c("Restrictions on CSOs (sum of net changes, 5yrs)"
                                ,"Restrictions on CSOs (level, 5yrs)"
                                , "Restrictions on CSOs, V-dem (sum of net changes, 5yrs)"
                                , "Restrictions on CSOs, V-dem (level, 5yrs)"
                                , "Conflict years (mean, 5yrs)"  
                                , "GDP (mean, 5yrs)"
                                , "Population (mean, 5yrs)")
          , type="text"
          , model.numbers=F
          , out="App_TableC4_Corruption_LevelRestrict_NewControls.txt"
          , title="Two-fixed effects linear models of corruption, different controls"
          , column.labels=c("Model C4a", "Model C4b", "Model C4c", "Model C4d")
          , dep.var.labels=c("Executive corruption index (V-Dem)")
          , notes = "(cluster-robust standard errors)")


#### Table C5 ####

# Standard Control Variables

set.seed(1234)

# DV: VDEM: Vote Buy and Clientelist Parties
formula1a <- as.formula( paste0("VDEM_votebuy ~ RESTRICT_COUNT_NOBAN_netl5 + RESTRICT_COUNT_NOBAN_nochange_netl5", controls,  " | cown  + year | 0 | cown") )
formula1b <- as.formula( paste0("VDEM_votebuy ~ RESTRICT_COUNT_NOBAN_nochange_netl5", controls,  " | cown  + year | 0 | cown") )
formula1c <- as.formula( paste0("VDEM_progrClientParties ~ RESTRICT_COUNT_NOBAN_netl5 + RESTRICT_COUNT_NOBAN_nochange_netl5", controls,  " | cown  + year | 0 | cown") )
formula1d <- as.formula( paste0("VDEM_progrClientParties ~ RESTRICT_COUNT_NOBAN_nochange_netl5", controls,  " | cown  + year | 0 | cown") )

# Model SG1
model1a <- felm(formula = formula1a,
                data = subset(data, year>=2000 & year<=2016 & BANNED==0)  )
summary(model1a)

model1b <- felm(formula = formula1b,
                data = subset(data, year>=2000 & year<=2016 & BANNED==0)  )
summary(model1b)

model1c <- felm(formula = formula1c,
                data = subset(data, year>=2000 & year<=2016 & BANNED==0)  )
summary(model1c)

model1d <- felm(formula = formula1d,
                data = subset(data, year>=2000 & year<=2016 & BANNED==0)  )
summary(model1d)

# Formula
formula1e <- as.formula( paste0("VDEM_votebuy ~ VDEM_CSOrepress_netl5 + VDEM_CSOrepress_nochange_netl5", controls,  " | cown  + year | 0 | cown + year") )
formula1f <- as.formula( paste0("VDEM_votebuy ~ VDEM_CSOrepress_nochange_netl5", controls,  " | cown  + year | 0 | cown + year") )
formula1g <- as.formula( paste0("VDEM_progrClientParties ~ VDEM_CSOrepress_netl5 + VDEM_CSOrepress_nochange_netl5", controls,  " | cown  + year | 0 | cown + year") )
formula1h <- as.formula( paste0("VDEM_progrClientParties ~ VDEM_CSOrepress_nochange_netl5", controls,  " | cown  + year | 0 | cown + year") )

# Model SG1
model1e <- felm(formula = formula1e,
                data = data)
summary(model1e)

model1f <- felm(formula = formula1f,
                data = data)
summary(model1f)

model1g <- felm(formula = formula1g,
                data = data)
summary(model1g)

model1h <- felm(formula = formula1h,
                data = data)
summary(model1h)


stargazer(model1a, model1b, model1c, model1d
          ,covariate.labels = c("Restrictions on CSOs (sum of net changes, 5yrs)"
                                ,"Restrictions on CSOs (level, 5yrs)"
                                , "Vertical accountability (mean, 5yrs)"
                                , "Horizontal accountability (mean, 5yrs)"
                                , "Conflict years (mean, 5yrs)" 
                                , "GDP (mean, 5yrs)"
                                , "Population (mean, 5yrs)")
          , type="text"
          , model.numbers=F
          , out="App_TableC5_Clientelism_LevelRestrict.txt"
          , title="Two-fixed effects linear models of clientelism"
          , column.labels=c("Model C5a", "Model C5b", "Model C5c", "Model C5d" )
          , dep.var.labels= c("Vote buying", "Clientelist Parties")
          , notes = "(cluster-robust standard errors)")

stargazer(model1e, model1f, model1g, model1h
          ,covariate.labels = c(  "Restrictions on CSOs, V-dem (sum of net changes, 5yrs)"
                                  , "Restrictions on CSOs, V-dem (level, 5yrs)"
                                  , "Vertical accountability (mean, 5yrs)"
                                  , "Horizontal accountability (mean, 5yrs)"
                                  , "Conflict years (mean, 5yrs)" 
                                  , "GDP (mean, 5yrs)"
                                  , "Population (mean, 5yrs)")
          , type="text"
          , model.numbers=F
          , out="App_TableC5cont_Clientelism_LevelRestrict.txt"
          , title="Two-fixed effects linear models of clientelism cont."
          , column.labels=c("Model C5e", "Model C5f", "Model C5g", "Model C5h")
          , dep.var.labels= c("Vote buying", "Clientelist Parties")
          , notes = "(cluster-robust standard errors)")


#### Table C6 ####

# Redefine Control Variables
controls_new <- paste0(" + Control_UCDP_ongoing_conflict_netl5 + Control_WB_controls_GDP_log_netl5 + Control_WB_controls_population_log_netl5")

set.seed(1234)

# DV: VDEM: Vote buying and clientelist parties
formula1a <- as.formula( paste0("VDEM_votebuy ~ RESTRICT_COUNT_NOBAN_netl5 + RESTRICT_COUNT_NOBAN_nochange_netl5", controls_new,  " | cown  + year | 0 | cown") )
formula1b <- as.formula( paste0("VDEM_votebuy ~ RESTRICT_COUNT_NOBAN_nochange_netl5", controls_new,  " | cown  + year | 0 | cown") )
formula1c <- as.formula( paste0("VDEM_progrClientParties ~ RESTRICT_COUNT_NOBAN_netl5 + RESTRICT_COUNT_NOBAN_nochange_netl5", controls_new,  " | cown  + year | 0 | cown") )
formula1d <- as.formula( paste0("VDEM_progrClientParties ~ RESTRICT_COUNT_NOBAN_nochange_netl5", controls_new,  " | cown  + year | 0 | cown") )

# Model SG1
model1a <- felm(formula = formula1a,
                data = subset(data, year>=2000 & year<=2016 & BANNED==0)  )
summary(model1a)

model1b <- felm(formula = formula1b,
                data = subset(data, year>=2000 & year<=2016 & BANNED==0)  )
summary(model1b)

model1c <- felm(formula = formula1c,
                data = subset(data, year>=2000 & year<=2016 & BANNED==0)  )
summary(model1c)

model1d <- felm(formula = formula1d,
                data = subset(data, year>=2000 & year<=2016 & BANNED==0)  )
summary(model1d)

stargazer(model1a, model1b,model1c, model1d
          ,covariate.labels = c("Restrictions on CSOs (sum of net changes, 5yrs)"
                                ,"Restrictions on CSOs (level, 5yrs)" 
                                , "Conflict years (mean, 5yrs)"  
                                , "GDP (mean, 5yrs)"
                                , "Population (mean, 5yrs)")
          , type="text"
          , model.numbers=F
          , out="App_TableC6_Clientelism_LevelRestrict_NewControls.txt"
          , title="Two-fixed effects linear models of clientelism, different controls"
          , column.labels=c("Model C6a", "Model C6b","Model C6c", "Model C6d")
          , dep.var.labels= c("Vote buying", "Clientelist Parties")
          , notes = "(cluster-robust standard errors)")

# Formula
formula1e <- as.formula( paste0("VDEM_votebuy ~ VDEM_CSOrepress_netl5 + VDEM_CSOrepress_nochange_netl5", controls_new,  " | cown  + year | 0 | cown + year") )
formula1f <- as.formula( paste0("VDEM_votebuy ~ VDEM_CSOrepress_nochange_netl5", controls_new,  " | cown  + year | 0 | cown + year") )
formula1g <- as.formula( paste0("VDEM_progrClientParties ~ VDEM_CSOrepress_netl5 + VDEM_CSOrepress_nochange_netl5", controls_new,  " | cown  + year | 0 | cown + year") )
formula1h <- as.formula( paste0("VDEM_progrClientParties ~ VDEM_CSOrepress_nochange_netl5", controls_new,  " | cown  + year | 0 | cown + year") )

# Model SG1
model1e <- felm(formula = formula1e,
                data = data)
summary(model1e)

model1f <- felm(formula = formula1f,
                data = data)
summary(model1f)

model1g <- felm(formula = formula1g,
                data = data)
summary(model1g)

model1h <- felm(formula = formula1h,
                data = data)
summary(model1h)

stargazer(model1e, model1f,model1g, model1h
          ,covariate.labels = c( "Restrictions on CSOs, V-dem (sum of net changes, 5yrs)"
                                 , "Restrictions on CSOs, V-dem (level, 5yrs)"
                                 , "Conflict years (mean, 5yrs)"  
                                 , "GDP (mean, 5yrs)"
                                 , "Population (mean, 5yrs)")
          , type="text"
          , model.numbers=F
          , out="App_TableC6cont_Clientelism_LevelRestrict_NewControls.txt"
          , title="Two-fixed effects linear models of clientelism, different controls (cont.)"
          , column.labels=c("Model C6e", "Model C6f","Model C6g", "Model C6h")
          , dep.var.labels= c("Vote buying", "Clientelist Parties")
          , notes = "(cluster-robust standard errors)")
